Mon. Not. R. Astron. Soc. 000, 1-11 (2009) Printed 5 January 2011 (MN style file v2.2) 



Models of our Galaxy II 

James Binney 1 * and Paul McMillan 1 

1 Rudolf Peierls Centre for Theoretical Physics, Keble Road, Oxford 0X1 3NP, UK 
Draft, September 8, 2010 



ABSTRACT 

Stars near the Sun oscillate both horizontally and vertically. In Paper I it was 
assumed that the coupling between these motions can be modelled by determining the 
horizontal motion without reference to the vertical motion, and recovering the coupling 
between the motions by assuming that the vertical action is adiabatically conserved 
as the star oscillates horizontally. Here we show that, although the assumption of 
adiabatic invariance works well, more accurate results can be obtained by taking the 
vertical action into account when calculating the horizontal motion. We use orbital tori 
to present a simple but fairly realistic model of the Galaxy's discs in which the motion 
of stars is handled rigorously, without decomposing it into horizontal and vertical 
components. We examine the ability of the adiabatic approximation to calculate the 
model's observables, and find that it performs perfectly in the plane, but errs slightly 
away from the plane. When the new correction to the adiabatic approximation is used, 
the density, mean-streaming velocity and velocity dispersions are in error by less than 
10 per cent for distances up to 2.5 kpc from the Sun. The torus-based model reveals 
that at locations above the plane the long axis of the velocity ellipsoid points almost 
to the Galactic centre, even though the model potential is significantly flattened. This 
result contradicts the widespread belief that the shape of the Galaxy's potential can be 
strongly constrained by the orientation of velocity ellipsoid near the Sun. An analysis 
of individual orbits reveals that in a general potential the orientation of the velocity 
ellipsoid depends on the structure of the model's distribution function as much as on 
its gravitational potential, contrary to what is the case for Stackel potentials. We argue 
that the adiabatic approximation will provide a valuable complement to torus-based 
models in the interpretation of current surveys of the Galaxy. 

Key words: The Galaxy: disc - The Galaxy: kinematics and dynamics - The Galaxy: 
structure - galaxies: kinematics and dynamics 



1 INTRODUCTION 

Study of the structure of the Milky Way Galaxy is a major 
theme of contemporary astronomy. Our Galaxy is typical of 
the galaxies that dominate the current cosmic star-formation 
rate, so understanding it is of more than parochial interest. 
We believe that most of its mass is contributed by elemen- 
tary particles that have yet to be directly detected, but we 
have only weak constraints on the spatial density and kine- 
matics of these particles - we urgently need stronger con- 
straints on them. The Cold-Dark-Matter (CDM) cosmogony 
provides a very persuasive picture of how a galaxy such as 
ours formed, and we need to know how accurately this the- 
ory predicts the structure of the Galaxy. 

In view of these considerations, large resources have 
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been invested over the last decade in massive surveys of the 
stellar content of the Galaxy. The rate at which data from 
this observational effort becomes available will increase at 
least through 2020. Models of the Galaxy will surely play 
a key role in extracting science from the data, because the 
Galaxy is a very complex object and every survey is sub- 
ject to powerful observational biases. Consequently it is ex- 
tremely hard to proceed in a model-independent way from 
observational data to physical understanding. We are likely 
to achieve the desired physical understanding by comparing 
observational data with the predictions of models. 

It is useful to distinguish between kinematic and dy- 
namic models. A kinematic model specifies the spatial den- 
sity of stars and their kinematics at each point without 
asking whether a gravitational potential exists in which 
the given density distribution and kinematics constitute a 
steady state. Bahcall & Soneira (1984) pioneered kinematic 
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models, and recent versions include Galaxia (Sharma et al. 
2010). The science of constructing dynamical models is still 
in its infancy. The Besangon model (Robin et al. 2003) has 
a dynamical element to it in that in it the disc's density pro- 
file perpendicular to the plane is dynamically consistent with 
the corresponding run of dispersion of velocities perpendic- 
ular to the plane. Schonrich & Binney (2009) and Binney 
(2010) (hereafter Paper I) offer models that are more thor- 
oughly dynamical. These models adopt a plausible model of 
the Galaxy's gravitational potential $(-R, z), in which there 
are substantial contributions to the local acceleration from 
a disc, a bulge and a dark halo, all assumed to be axisym- 
metric. They assume that motion parallel to the plane is 
to some degree decoupled from motion perpendicular to the 
plane. Specifically, the vertical motion is governed by the 
time-dependent potential 

*,(*;«) = S[Ji(t),*] (1) 

where R(t) is the radius at time t that one obtains by 
assuming that the radial motion is governed by the one- 
dimensional effective potential 

* R (R) = *(R,0) + §- (2) 

Paper I assumed that the time-dependence of the potential 
(1) is slow enough that the action J z of vertical motion is 
constant. It justified this assumption by referring to Fig- 
ure 3.34 in Binney & Tremaine (2008) (hereafter BT08), 
which shows that the boundaries of one particular orbit are 
fairly well recovered by the adiabatic approximation (here- 
after AA). In this paper we explore the validity of the AA 
much more extensively. Our other goal is to present a model 
of the Galactic disc that is not reliant on the AA. This model 
dispenses with the assumption that the R and z motions are 
decoupled by using numerically synthesised orbital tori. 

The paper is organised as follows. In Section 2 the va- 
lidity of the AA is tested on typical orbits. In Section 3 we 
explain the general principles of torus modelling and why 
we believe this technique will prove a valuable tool for the 
interpretation of observational data. We define the torus- 
based model of the Galactic discs that we will use to test 
the accuracy of observables obtained from the AA, and we 
summarise the methods used to extract observables when 
the model is based on (i) tori, and (ii) the adiabatic approx- 
imation. In Section 4 we compare the model's observables 
with estimates of them obtained from the AA. In Section 5 
we examine the tilt of the velocity ellipsoid near the Sun, 
which cannot be computed from the AA, and show that its 
long axis points towards the Galactic centre even though 
the potential is significantly flattened. Section 6 sums up 
and looks to the future. An appendix explains how some 
important Jacobians are calculated for a torus model. 



2 VALIDITY OF THE ADIABATIC 
APPROXIMATION 

Each panel of Fig. 1 shows an orbit in the meridional plane 
in the gravitational potential of a Miyamoto-Nagai Galaxy 
model (Miyamoto & Nagai 1975) with scale-length ratio 
b/a = 0.2: Figure 2.7 of BT08 shows that this model has 
a prominent disc. From the Fourier transforms of the time 
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Figure 1. Two orbits in the meridional plane of a Miyamoto- 
Nagai model with scale-length ratio b/a = 0.2. In both cases the 
angular momentum about the symmetry axis is L z = y/GM/a, 
where M is the model's mass. In units of V GMa, the actions 
of the upper and lower orbits are (J r ,Jz) = (0.109,0.067) and 
(0.127,0.022). The numbers above each panel are the values of 
J z of obtained by following the motion of particles dropped from 
the upper left and upper right corners of the orbit in the one- 
dimensional potential ^ z (z) = &(R, z) — &(R, 0), and the corre- 
sponding vertical energies. The nearly straight full lines show the 
AA to the orbit when J z is set to the average of the values at top 
left and the radial action takes its true value. The dashed lines 
show the boundary yielded by the AA when L z is replaced in the 
effective potential by L z + J z . 

series R(t) z(t) on these orbits (Binney & Spergel 1984) we 
find that in units of V GMa their actions are (J r ,Jz) = 
(0.109,0.067) and (0.127,0.022), respectively. One can also 
estimate the vertical actions J z of these orbits by dropping 
particles from points on their upper edges, and following 
their motion in the one-dimensional potential (1) with R 
frozen at its current value - see equation (9) below. The 
numbers at top of each panel show the values obtained for 
J z in the same units when particles are dropped from the 
top left and top right corners of the orbit; the values are 
displaced from the true value by < 7%. The corresponding 
vertical energies E z are also shown at the top of each panel; 
they differ by more than a factor 2. Thus the AA does provide 
a fairly good guide to how the vertical motion is influenced 
by the radial motion. 

The nearly straight full lines in Fig. 1 show the outlines 
of the approximate orbits that we obtain from the AA by set- 
ting J r to its true value and J z to the average of the values 
given at the top of the panel. The shape of each approximate 
orbit is reasonable, although the left and right edges are 
straight rather than curved, but the orbit is clearly displaced 
to smaller radii with respect to the true orbit. This differ- 
ence reflects the fact that the vertical motion contributes to 
the centrifugal potential alongside the azimuthal motion. By 
assuming that the radial motion is governed by the effective 
potential (2) in which L z occurs rather than the total angu- 
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Figure 2. Surfaces of section R/a = 1.9, 2.4 and 3 for the orbit shown in the top panel of Fig. 1. The points are for the numerically 
integrated orbit, while the curves are obtained from the AA. 



lar momentum L, we have under-estimated the centrifugal 
potential. Consequently, we predict that the orbit lives at 
smaller radii than it really does. 

In a spherical potential, the total angular momentum 
is related to L z and J z by L = \L Z \ + J z (e.g. BT08 §3.5.2) 
and the radial motion is governed by an effective potential 
in which the centrifugal component is L 2 /2r 2 , where r 2 = 
R 2 + z 2 . Consequently, an obvious strategy for improving 
the predictions of the AA is to replace L z in the effective 
potential by L + J z . The dashed lines in Fig. 1 show the 
effect of replacing L z by 

C z = \L Z \ + 1 J Z , (3) 

with 7 = 1. Both orbits are now quite closely modelled. 

If we calculate the radial action of a given phase-space 
point (x, v) using L z rather than C z in the effective poten- 
tial, the value we obtain is too large when (x, v) lies near 
apocentre, because the star moves in an effective potential 
that has its minimum at a radius that is too small. Con- 
versely, when (x, v) lies near pericentre, our estimate of J r 
is too small if we use L z . Since the df decreases with increas- 
ing J r , the use of the less accurate effective potential leads 
to phase-space points near pericentre being over-weighted 
relative to points near apocentre, and this in consequence 
shifts the predicted distribution in to large values. Hence, 
replacing L z in the effective potential with C z for suitably 
chosen 7 can usefully improve the accuracy of results ob- 
tained with the AA. 

The points in Fig. 2 show the consequents of the up- 
per orbit of Fig. 1 in three surfaces of section that are ob- 
tained by noting z and p z when the star crosses the line 
R — constant in the meridional plane with pr > 0. The 
curves in each panel show the dependence of p z on z along 
the one-dimensional orbit in the potential <fr(R, z) with R 
fixed at the appropriate value when the action J z is set to 
the average of the values given above the top panel of Fig. 1. 
The agreement between the curves yielded by the AA and the 
numerical consequents is on the whole good. In the left and 
central panels we see that while the curves have reflection 
symmetry in p z = the consequents do not. This is because 
the surface of section is for pn > 0, and when the star is 
moving outwards, it is likely to be moving upwards when 
z > and downwards when z < 0. As we will discuss in 
Section 5, when a galaxy is formed out of such orbits, this 
z-dependent correlation between pn and p z causes the prin- 
cipal axes of the velocity ellipsoid to become inclined to the 
R, z axes at \z\ > 0. The AA is unable to capture this aspect 



of the dynamics and will always yield a velocity ellipsoid 
that is aligned with the R, z axes. 

The panel on the extreme right shows that at large radii 
the AA underestimates the maximum height z max reached by 
a star, although at most values of z it predicts p z with good 
accuracy. 

The analogue of Fig. 2 for the orbit shown in the lower 
panel of Fig. 1 shows smaller offsets between the numerical 
consequents and the predictions of the AA, because the latter 
works best for small vertical amplitudes. 

With vt denoting the tangential speed, the centrifugal 
potential is v 2 /r 2 . In a separable potential, the time aver- 
age of the ith component of velocity is related to the ith 
frequency and action by (v 2 ) = QiJi, so when we replace v 2 
by its time average the centrifugal potential becomes 

\L z \ + (n z /^)j z 

WTz-z • (4) 

The standard AA underestimates the centrifugal potential 
by neglecting the term proportional to J z in the numera- 
tor, and partially compensates by neglecting the z 2 in the 
denominator. This neglect of z 2 must be responsible for the 
fact that we find the optimum value of 7 to be unity rather 
than 2, which is a typical value of Q. z /Q.$, for disc stars at 
Ro in plausible Galactic potentials. However, in an unreal- 
istically flat potential, larger values of 7 prove optimal. For 
example, when the potential is that of a razor-thin expo- 
nential disc and there is no contribution from a bulge or a 
dark halo, we find 7 ~ 1.9 is required. Even in this case 7 is 
smaller than the typical value of Qz/fi^, on account of the 
neglect of z in the denominator. 



3 A MODEL BASED ON ORBITAL TORI 

The classical approach to modelling globular clusters starts 
by positing an analytic form for the distribution function 
(df) and then calculating the density distribution and kine- 
matics that are implied by this df. Thus globular clusters 
have been successfully modelled with dfs of the King-Michie 
form (e.g. BT08 §4.3). This approach can be extended to disc 
galaxies. For example Rowley (1988) modelled SO galaxies 
with distribution functions of the form f(E,L z ), where E 
and L z are, respectively, orbital energy per unit mass and 
and angular momentum per unit mass about the symme- 
try axis. Unfortunately, such a simple distribution function 
cannot successfully model the Galaxy, because it predicts 
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equal velocity dispersions an. and o z in the radial and verti- 
cal directions, while observations show that ob, — 1.1a z (e.g. 
Aumer & Binney 2009). A successful df for the Galaxy must 
depend, explicitly or otherwise, on the vertical action J z . 

Given that the df will depend on two of an orbit's three 
actions J r , J z and L z , substantial advantages arise from em- 
ploying J r as the other argument of the df in place of E. 
For this reason Paper I studied Galaxy models in which each 
stellar component had a df that was an analytic function 
of the three actions. It used the AA to calculate observables 
from the df. The purpose of this section is to compare ob- 
servables obtained in this way to those obtained without 
invoking the A A but instead using orbital tori. 

3.1 General principles of torus modelling 

Orbital tori are the three-dimensional surfaces in six- 
dimensional phase space on which individual orbits move. 
They are the building blocks from which Jeans' theorem as- 
sures us that any equilibrium model can be built. Each torus 
is characterised by a set of three actions J = (J r , L z , J z ) and 
therefore corresponds to a point in action space. We build a 
galaxy model by assigning a weight to each torus. 

We obtain tori as the images of analytic tori under a 
canonical transformation. The tori used here are defined 
by the angle-action coordinates of the isochrone potential 
(e.g. BT08 §3.5). Given actions J, the computer constructs a 
canonical transformation that maps the analytic torus with 
actions J into the required torus by adjusting the coefficients 
in a trial generating function so as to minimise the rms vari- 
ation of the Galactic Hamiltonian on the image torus. Once 
this has been done, we have analytic expressions for the 
phase-space coordinates [x(0), v(0)] as functions of the an- 
gle variables which control the orbital phase. On a given 
torus, the phase-space coordinates (x, v) are 27r-periodic 
functions of each The torus-fitting program also returns 
the values of the torus's characteristic frequencies Q», so we 
can determine the motion of a star using 8(t) — 0(0) + fit. 
For a fuller summary of how orbital tori are constructed 
and references to the papers in which torus dynamics was 
developed, see McMillan & Binney (2008). 

Torus modelling is best understood as an extension of 
Schwarzschild modelling (Schwarzschild 1979), which has 
been successfully used in many studies of the dynamics of 
external galaxies (e.g. Gebhardt et al. 2003; Krajnovic et 
al. 2005). A Schwarzschild model is constructed by assign- 
ing weights to each orbit in a "library" of orbits. The orbit 
library is assembled by integrating the equations of motion 
in the given potential for a sufficient time, and noting the 
fraction of its time that the orbiting particle spends in each 
bin in the space of observables. Then a non-negative weight 
Wi is chosen for each orbit such that the data are consistent 
with the model's predictions. In torus modelling orbits are 
replaced by tori, which are essentially equivalence classes 
of orbits that differ from one another only in phase, and a 
Runge-Kutta integrator is replaced by the torus-fitting code. 
Whereas orbits are defined by their six-dimensional initial 
conditions, tori are defined by their actions J. 

Replacing numerically integrated orbits with orbital 
tori brings the following advantages 

(i) The phase-space density of orbits becomes known be- 



cause tori have prescribed actions and the six-dimensional 
phase-space volume occupied by orbits with actions in d 3 J 
is t — (27r) 3 d 3 J. Knowledge of the phase-space density of 
orbits allows one to convert between orbital weights and the 
value taken by the df on an orbit. 

(ii) On account of the above result, there is a clean and 
unambiguous procedure for sampling orbit space and relat- 
ing the weights of individual tori to the value that the df 
takes on them - see below. The choice of initial conditions 
from which to integrate orbits for a library is less straight- 
forward because the same orbit can be started from many 
initial conditions, and when the initial conditions are sys- 
tematically advanced through six-dimensional phase space, 
the resulting orbits are likely at some point to cease explor- 
ing a new region of orbit space and start resampling a part 
of orbit space that is already represented in the library. On 
account of this effect, it is hard to relate the weight of an 
orbit to the value taken by the df on it (but see Hafher et 
al. 2000; Thomas et al. 2005). 

(iii) There is a simple relationship between the distribu- 
tion of stars in action space and the observable structure and 
kinematics of the model; as explained in §4.6 of BT08, the 
observable properties of a model change in a readily under- 
stood way when stars are moved within action space. The 
simple relationship between the observables and the distri- 
bution of stars in action space enables us to infer from the 
observables the form of the df /(J), which is nothing but 
the density of stars in action space. 

(iv) From a torus one can readily find the velocities that 
a star on a given orbit can have when it reaches a given 
spatial point x. By contrast a finite time series of an orbit 
is unlikely to exactly reach x, and searching for the time at 
which the series comes closest to x is laborious. Moreover, 
several velocities are usually possible at a given location, and 
a representative point of closest approach must be found for 
each possible velocity. 

(v) An orbital torus is represented by of order 100 num- 
bers while a numerically-integrated orbit is represented ei- 
ther by some thousands of six-dimensional phase-space lo- 
cations, or by a similar number of occupation probabilities 
within a phase-space grid. 

(vi) The numbers that characterise a torus are smooth 
functions of the actions J. Consequently tori for actions that 
lie between the points of any action-space grid can be con- 
structed by interpolation on the grid. Interpolation between 
time series is not practicable. 

(vii) Schwarzschild and torus models are zeroth-order, 
time-independent models which differ from real galaxies 
by suppressing time-dependent structure, such as ripples 
around early-type galaxies (Malin & Carter 1980; Quinn 
1984; Schweizer & Seitzer 1992), and spiral structure or 
warps in disc galaxies. Since the starting point for pertur- 
bation theory is action-angle variables (e.g. Kalnajs 1977), 
in the case of a torus model one is well placed to add time- 
dependent structure as a perturbation. Kaasalainen (1995) 
showed that classical perturbation theory works extremely 
well when applied to torus models because the integrable 
Hamiltonian that one is perturbing is typically much closer 
to the true Hamiltonian than in classical applications of per- 
turbation theory (Gerhard & Saha 1991; Dehnen & Gerhard 
1993; Weinberg 1994), in which the unperturbed Hamilto- 



Models of our Galaxy II 5 



Table 1. Parameters of the DF. 

Disc JJd/kpc cr r o/kms _1 o^o/kms -1 Lo/kpc km s — 1 



Thin 
Thick 



2.4 
2.5 



27 
48 



20 
44 



10 
10 



nian arises from a potential that is separable (it is generally 
either spherical or plane-parallel). 



3.2 Choice of the DF 

For the comparison of results obtained with and without 
the AA it is appropriate to study a model that has a very 
simple DF. Specifically we represent both the thin and the 
thick discs with a DF that is quasi-isothermal in the sense of 
Paper I: 



/(J r ,L z , J z ) = U r (Jr,L z ) X 



27rcr| 



where 



fa r (Jr, L z ) 



[1 + tanh(L z /L )]e~ 



(5) 



(6) 



Here Q(L Z ) is the circular frequency for angular momen- 
tum L z , k(L z ) is the radial epicycle frequency and v(L z ) 
is its vertical counterpart. £(L Z ) = E e~^ R ~ Rc ^ /Rd is the 
(approximate) radial surface-density profile, where R C {L Z ) 
is the radius of the circular orbit with angular momentum 
L z . The factor l+tanh(L z /Lo) in equation (6) is there to ef- 
fectively eliminate stars on counter-rotating orbits and the 
value of Lq is unimportant provided it is small compared 
to the angular momentum of the Sun. In equations (5) and 
(6) the functions a z (L z ) and ay(L z ) control the vertical and 
radial velocity dispersions. The observed insensitivity to ra- 
dius of the scaleheights of extragalactic discs motivates the 
choices 

ar(L z ) = * r0 ei (R °- R ^/ R « 

* z (L z ) = a z0 e q< - Ro - R c )/R «, (7) 

where q = 0.45 and cr r o and a z o are approximately equal to 
the radial and vertical velocity dispersions at the Sun. We 
take the DF of the entire disc to be the sum of a DF of the 
form (5) for the thin disc, and a similar DF for the thick 
disc, the normalisations being chosen so that at the Sun the 
surface density of thick-disc stars is 23 per cent of the total 
stellar surface density. Table 1 lists the parameters of each 
component of the DF. 

There are two main differences between the DF we use 
here and that used in Paper I: (i) Paper I used the actual ver- 
tical frequency Q Z (J) in equation (5) while here we use the 
vertical epicycle frequency v{L z ). This substitution is nec- 
essary because for large J, Q. z tends to zero so fast that the 
product Q z J z can decrease as J z — > oo, leading to unphys- 
ical results when Q Z J Z appears in the DF as the argument 
of an exponential, (ii) In the interests of simplicity the thin 
disc is here represented by a single quasi-isothermal compo- 
nent whereas in Paper I it was represented it by a sum of 
quasi-isothermals, one for stars of each age. 

Any serious attempt to fit a real stellar catalogue must 
distinguish between stars of different ages, and different 
metallicities, because the colours and luminosities of stars 



are very much functions of age and metallicity, so the 
chances of a star entering a catalogue depend on its age 
and metallicity. Consequently, by lumping together all thin- 
disc stars regardless of age we forgo the opportunity to fit a 
real stellar catalogue in a detailed way. Nonetheless, we shall 
require our DF to reproduce an observational density profile 
to demonstrate that even our unrealistically simple DF has 
sufficient flexibility to reproduce given data to reasonable 
precision. 

The physical properties of the model are jointly deter- 
mined by the DF and the gravitational potential <E>(i?, z). 
Ultimately it will be necessary to require that the Galaxy's 
DF be consistent with <E> in the sense that the density of 
matter that the DF predicts generates However, before 
the question of dynamical self-consistency can be addressed, 
one must not only specify the DF of dark matter (which is 
believed to contribute about half the gravitational force on 
the Sun) but also distinguish carefully between the masses 
of stars and their luminosities in the wavebands in which 
they are observed. In practice the latter can be done only if 
one has specified the Galaxy's star-formation and metal- 
enrichment history. This enterprise goes far beyond the 
scope of the present paper; it will be addressed in subse- 
quent papers in this series, which will explain the impor- 
tance of comparing models to data in the space of observ- 
ables, such as apparent magnitudes, parallaxes and proper 
motions, rather than the space of physical variables such 
as (x, v) used here. The purpose of this paper is merely to 
lay the foundations for such an exercise, which we expect 
to give first insights into the DF of dark matter. Here we 
take the view that our DF weights stars by their luminosity 
rather than their mass, and assume that <3? is the potential of 
Model 2 of Dehnen & Binney (1998) modified to have thin- 
and thick-disc scaleheights of 360 pc and 1 kpc (Table 2) . In 
this model the disc contributes 60 per cent of the gravita- 
tional force on the Sun, with dark matter contributing most 
of the remaining force. 



3.3 Modelling procedures 

Together the DF and the potential specify the probability 
density of stars in phase space. The simplest way to derive 
the model's physical characteristics from this probability 
density is to obtain a discrete realisation of the probabil- 
ity density by Monte-Carlo sampling. The model's physical 
characteristics are then obtained by binning the realisation's 
stars. The DF specified by equations (5) and (6) can be an- 
alytically integrated over J z and J r to obtain the marginal 
distribution in L z , so we can obtain a discrete realisation 
of this DF by successively sampling one-dimensional pdfs in 
L z , J r and J z . The results presented below are typically 
obtained with ~ 200 000 tori. 

Once we have a torus library, a discrete realisation of 
the Galaxy is obtained by repeatedly choosing a torus from 
the library at random, then choosing each angle variable 
uniformly within (0, 2ir), and using the functions returned 
by the torus-generating software to to determine (x, v) from 
the given values of J and 8. 

When the AA is used, model construction proceeds 
rather differently: then given (x, v) one determines J r and 
J z by the following steps. 
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Table 2. Parameters of the potential 
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Figure 3. Top: density as a function of distance from the mid- 
plane at the solar radius. Bottom: mean streaming velocity as 
a function of distance from the midplane. The full black curves 
show the predictions of the full torus model; the blue curves are 
obtained from the A A with 7 = 0; the dotted red curves show the 
effect of setting 7 = 1. 



• Evaluate the vertical and radial energies 



E z = \v 2 z + <S> z {z) 



Er = \v R + *r(R), 

where ^ z and ^r are defined by equations (1) and (2). 
• Evaluate the actions from 



(8) 



Jz = - dzv z (E z , z) 

TV ' 



dRv R (E R ,R), 



(9) 



where z max is defined by ^/ z (z inax ) = E z and R a and R p are 
defined by = E R . 

These steps make it straightforward to evaluate the df at 
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Figure 4. The distributions for local stars of the radial, azimuthal 
and vertical components of velocity after marginalising over the 
other two components. The full black curves are for the torus- 
based model while the broken curves are obtained using the AA 
with 7 = (dashed blue) and 7 = 1 (dotted red). The curves 
overlie one another too closely to be clearly distinguishable. 



an arbitrary point (x, v), and thus derive the model's phys- 
ical properties by numerically integrating the df, times any 
power of Vi, over velocity space. The quantities such as stel- 
lar number density ^(x) and velocity dispersion (u z ) 1//2 (x) 
are then continuous functions of their arguments. In the 
absence of the AA, an iterative procedure such as that de- 
scribed by McMillan & Binney (2008) is required to deter- 
mine J(x, v), and the torus-modeeling procedure avoids this 
procedure by choosing J not (x, v). The price we pay for 
starting with J is discreteness, and the necessity of estimat- 
ing f(x), etc, by binning stars. 
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4 COMPARISONS 

Fig. 3 shows the density of stars (upper panel) and the mean- 
streaming velocity (lower panel) as functions of \z\ at the so- 
lar radius. The full curves show results from the torus model, 
while the dashed and dotted curves show results obtained 
from the AA using two values of the parameter 7 defined 
by eq. (3): 7 = (dotted red) and 7=1 (dashed blue). 
Also shown in the upper panel are the seminal data points 
of Gilmore & Reid (1983), which led to the identification of 
the thick disc. Since our df provides a reasonable fit to these 
points, it may be close to the actual df for turnoff stars. The 
AA recovers the density profile of the torus model to good 
accuracy for either value of 7. 

The lower panel in Fig. 3 shows how the mean stream- 
ing speed (v,f,) is predicted to fall with \z\. The agreement 
between the torus model and the model based on the AA 
with 7=1 (dotted red) is excellent for \z\ < 1.6 kpc but at 
larger heights the AA has a systematic tendency to under- 



estimate {v$). On account of the problem discussed in Sec- 
tion 2 apropos Fig. 1, the AA with 7 = over-estimates the 
mean-streaming speed by a few kms -1 for \z\ >0.7kpc. 

Fig. 4 shows the distributions of the radial ([/), tangen- 
tial (V) and vertical (W) components of velocity in a small 
volume akin to the solar neighbourhood. Since the full black 
curves from the torus model coincide with the dashed blue 
curves from the AA with 7 = 0, the AA reproduces the torus 
results to high precision. The results obtained on setting 
7 = 1 are plotted as a dotted red curve but overlie the other 
curves too closely to be clearly distinguishable. Fig. 5 shows 
that the AA is equally successful in recovering the distribu- 
tions of U, V and W at a smaller Galactocentric radius, 
6.5 kpc. The insensitivity to 7 of the velocity distribution in 
the plane arises because these distributions are dominated 
by rather nearly circular orbits, and for these orbits J z is so 
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Figure 7. The radial, azimuthal and vertical velocity dispersions 
as functions of distance from the midplane. The full black curve 
shows results from the torus model, while the dashed blue line is 
obtained using the AA with 7 = 0; the dotted red curves show the 
effect of setting 7 = 1. 

much smaller than L z that adding J z to L z barely changes 
the numerical value. 

Fig. 6 shows the U, V and W distributions at R = 8 kpc 
and 1.5 kpc away from the midplane. Systematic differences 
between the predictions of the AA and the results of the full 
torus model are now evident. For either value of 7, the AA 
yields a distribution of U that is slightly too broad. When 
7 = 0, the AA gives a distribution in W (dashed blue curve) 
that is too sharply peaked, but this fault is nicely corrected 
by setting 7 = 1 (dotted red curve) because increasing 7 
moves orbits of given L z outwards, and by virtue of equation 
(7), the smaller an orbit's value of L z , the faster it is likely 
to move vertically. As expected, with 7 = the AA yields 
a distribution in V that is offset from that of the full torus 
model by ~ 8kms _1 towards higher velocities (dashed blue 
curve). This offset is largely cured by setting 7=1 (dotted 
red curve). 




7 8 9 10 11 

R 



Figure 8. Three effective potentials for the orbit in the Galaxy 
model that has actions (J r ,L z ,J z ) = (0.05,1,0.05) times the 
angular momentum of the circular orbit at Rq . This orbit extends 
to 3 kpc above the plane. The dotted and dashed curves show the 
naive effective potential (2) with L z replaced by \L Z \ + -yj z and 
7 = in the dotted case and 7 = 1 in the dashed case. The full 
curve shows the effective potential derived from a time-average of 
the radial force. 



Fig. 7 shows the variation with \z\ of the radial, tangen- 
tial and vertical velocity dispersions. The two planar veloc- 
ity dispersions are accurately reproduced for |z|<1.3kpc. 
At greater heights the AA over-estimates the dispersions, 
most strikingly so in the case of a^. The excessive value of 
<j$ is clearly associated with the tendency of the red curve 
in the middle panel of Fig. 6 to lie above the black one at 
i)0 ~ 50kms -1 . As the top and middle panels of Fig. 6 
suggest, (ujj) 1//2 and {(v^ — (^)) 2 ) 1//2 are at any height re- 
markably insensitive to the value of 7. 

From the bottom panel of Fig. 6 we anticipate that 
increasing 7 will increase {v z } at large \z\ and indeed the 
bottom panel of Fig. 7 shows that increasing 7 from zero to 
unity increases a z at all z, but particularly at large z. This 
result arises because increasing 7 shifts orbits with large J z 
outwards, and since /(J) is a strongly decreasing function 
of I J I, this outwards shift raises the density of stars with 
large J z at a given location. At |z|>lkpc setting 7=1 
increases the accuracy of a z , while at smaller values of \z\, 
a slight deterioration in accuracy results. The bottom panel 
of Fig. 6 hints that the full curve in Fig. 7 may lie slightly 
too low as a result of poor sampling by the torus model of 
orbits with large \v z \, so the results from the AA with 7=1 
may be more accurate than appears to be the case. 

The dotted and dashed curves in Fig. 8 show the ef- 
fective potential (2) for a typical solar-neighbourhood orbit 
when 7 = and 7 = 1, respectively. The full curve shows an 
effective potential for this orbit that was obtained by first 
evaluating the time-average of 9$ /dR- L 2 Z /R 3 at each value 
of R visited by the orbit, and then integrating the resulting 
function of R. We see that with L z replaced by \L Z \ + J Z the 
simple effective potential (2) provides a good fit to the ef- 
fective potential obtained by time-averaging the radial force 
as a function of R. 
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arctan(z/R) 

Figure 9. The variation of the angle of tilt of the velocity el- 
lipsoid with respect to the plane versus the angle arctan(|z|/i?) 
between the plane and the line of sight from (R, z) to the Galac- 
tic centre. Since the full curve tracks the dotted line, the long 
velocity ellipsoid points almost straight at the Galactic centre. 




Figure 10. The distribution of the ratio u z /ur that gives the di- 
rection of a principal axis of the contribution of each torus to the 
velocity ellipsoid at the point (R, z) = (Ro, 1.5 kpc). The full his- 
togram weights each torus equally, while the dotted one weights 
them in proportion to their contributions to the density at the 
given point. Even the latter distribution is wide, so the orien- 
tation of the velocity ellipsoid depends quite sensitively on the 
weights assigned to orbits by the DF. 



5 TILT OF THE VELOCITY ELLIPSOID 

A popular diagnostic of the Galaxy's gravitational potential 
is the way in which the principal axes of the velocity ellipsoid 
tilt as one moves away from the plane (Siebert et al. 2008; 
Smith et al. 2009). As stated above, this phenomenon lies 
beyond the scope of the AA, but it can be determined from 
the torus model. Fig. 9 shows that this angle is nearly equal 
to the angle arctan(|«|/i?) between the plane and the line 
of sight from (R, z) to the Galactic centre. That is, in the 
vicinity of the Sun, the longest axis of the velocity ellipsoid 
is almost parallel to the radial vector r. 

This behaviour is that expected in a spherical potential, 
and Smith et al. (2009) argue that alignment of the velocity 
ellipsoid with spherical coordinates implies that the poten- 
tial is spherically symmetric. Our assumed gravitational po- 



tential is far from spherical because roughly half the radial 
force at the Sun is contributed by the disc, and the dark 
halo is itself flattened (axis ratio 0.8). The aspherical nature 
of the potential is reflected in the fact that the distribution 
of frequency ratios O z /f20 of the model's tori, has a me- 
dian close to 2, while in a spherical potential this ratio is 
inevitably unity. 

Although our model is not strictly speaking a counter 
example to the assertion of Smith et al. (2009) because 
we have established only that the velocity ellipsoid is ap- 
proximately aligned with spherical coordinates in the region 
around the Sun, it does suggest that one should examine 
more closely the reasoning in that paper. 

A key step in its argument is the assertion that if the 
DF is an even function of then the isolating integrals that 
constrain individual orbits are also. If the isolating integrals 
have this property, then the velocity ellipsoid provided by 
any DF will be aligned with spherical coordinates. In partic- 
ular the DF /(J) = 5(3 — Jo) corresponding to a single orbit 
will be radially aligned, so the matrix (viVj) x will be diag- 
onal in spherical coordinates, where (.) x implies the time 
average over instants when the star lies in some small vol- 
ume around x. 

As Fig. 1 illustrates, for stars on standard orbits in the 
meridional plane, four velocities are possible at a given point. 
Let these velocities be ±vi and ±V2. By time-reversal sym- 
metry, the probability that vi occurs is inevitably equal to 
the probability that — vi occurs, and similarly for ±V2. How- 
ever, it turns out that ±vi may occur more or less often 
than ±V2. In fact the probability of occurrence of is pro- 
portional to the Jacobian d(0)/d(x), which is a non-trivial 
function of 0, but one that can be determined for a numeri- 
cally constructed torus (see Appendix A). The four possible 
velocities at x correspond to the four values of 9 that bring 
the star to the given point x, and the values p\ and p2 taken 
by d(0)/d(x.) at these values of 8 depend on whether v is 
±vi or ±V2. Clearly we have 



(ViVj) x 



PlVliVlj +P2V2iV2j 
Pi +P2 



(10) 



Let u be the eigenvector of this matrix that lies closest to the 
r direction. Fig. 10 shows for the point (R, z) = (8, 1.5) kpc 
the distribution of the ratio u z /u R . The full histogram shows 
the distribution when each torus is given equal weight, while 
the dotted histogram shows the distribution when tori are 
correctly weighted by the contributions that they make to 
the stellar density at the given point. (Although the density 
of sampling in action space ensures that all tori make equal 
contributions to the stellar mass of the entire Galaxy, every 
torus has its own way of spreading its mass in space.) We see 
that the distribution of orientations of individual contribu- 
tions to the velocity ellipsoid is quite broad, even when the 
tori are correctly weighted. An examination of the depen- 
dence of u z /ua on J z reveals that orbits with larger values 
of J z make contributions that are aligned with r, while it 
is orbits with small J z that sometimes make contributions 
that are aligned nearly parallel to the plane. Since the DF 
specifies the relative weight of these variously oriented con- 
tributions, it controls the orientation of the final velocity 
ellipsoid at least as much as does the gravitational poten- 
tial. Consequently, only limited inferences about the nature 
of the potential can be drawn from observations of the ve- 
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locity ellipsoid yielded by the df that the Galaxy happens 
to have. 

The widespread belief that the shape of the potential 
can be inferred from the orientation of the velocity ellipsoid 
probably arises from studies of models that have Stackel po- 
tentials. For these potentials the Hamilton-Jacobi equation 
separates in an appropriate coordinate system (u,v), and 
the canonically conjugate momenta are functions of one co- 
ordinate only: p u {u), p v (v) (e.g. BT08 §3.5.3). Consequently, 
the coordinate directions are bisectors of the angles between 
vi and V2. Moreover, it turns out that for these potentials, 
d(9)/d(x) is the same for all four values of 8 so the co- 
ordinate directions are the eigenvectors of {viVj) x for every 
orbit that reaches x. Consequently, the velocity ellipsoid has 
to be oriented with the coordinate directions regardless of 
how orbits are weighted. 1 In a general potential, there is 
no universal coordinate system that describes the alignment 
of the egenvectors of {viVj) x , and the orientation of the fi- 
nal velocity ellipsoid very much depends on how orbits are 
weighted. 



6 CONCLUSIONS 

Dynamical models of the Milky Way will play an key role 
in the scientific exploitation of data from large surveys that 
are currently being undertaken. Models that are based on 
Jeans' theorem should be the most powerful tools for ex- 
tracting science from data, and amongst such models those 
that express the df as a function of the actions enjoy some 
very important advantages. 

The major obstacle to the use of Jeans' theorem in the 
context of the Galaxy is the lack of analytic expressions 
for three independent isolating integrals. Paper I presented 
models that use approximate expressions for the actions that 
rely on the adiabatic invariance of the vertical action J z . In 
Section 2 we tested the validity of this adiabatic approxima- 
tion (aa) by numerically integrating typical orbits. We found 
that the orbits' vertical dynamics is reproduced by the AA to 
remarkably good accuracy, but the motion in the plane is less 
accurately recovered because the naive AA under-estimates 
the strength of the centrifugal potential. This defect leads 
to the radial action derived for a given phase-space point 
(x, v) being over-estimated when the point lies near apocen- 
tre, and under-estimated when it lies near pericentre. Since 
t>0 is small at apocentre and large near pericentre and the 
df is a declining function of J r , the defect leads to the mean- 
streaming velocity {v$) being over-estimated. The problem 
can be largely resolved by replacing L z in the centrifugal 
potential by \L Z \ + ^J z with 7 a number of order unity. 

The more strongly flattened the potential is, the more 
accurate the AA becomes and the larger the value of 7 needs 
to be. For example, in the extreme case of vanishing dark 
halo 7 = 1.9 works well. 

In Section 3 we explained how to build a model Galaxy 
using orbital tori. Torus modelling is best considered an ex- 
tension of Schwarzschild modelling, which has long been 
a standard tool for the interpretation of data for external 

1 This conclusion can be obtained more simply by observing that 
the isolating integrals E and I3 upon which the DF must depend, 
are even functions of p u and p v . 



galaxies, both in connection with searches for massive black 
holes and attempts to understand how early-type galaxies 
were assembled. Torus modelling is a more powerful tech- 
nique principally because (i) it enables us to quantify or- 
bits by the values taken on them of essentially unique and 
physically easily understood isolating integrals, and (ii) it 
makes it easy to determine at what velocities a star will 
pass through any spatial point. We presented a df of ex- 
ceptional simplicity, which generates a reasonably realistic 
model of the Galaxy's discs. 

In Section 4 we examined in some detail observable 
quantities in this model when they are calculated from ei- 
ther the full torus machinery, or from the AA. We showed 
that in the plane, at both Ro and R = 6.5 kpc, the distribu- 
tions of all three components of velocity are reproduced to 
high accuracy by the AA, regardless of whether 7 is set to 
zero or unity. Away from the plane the velocity distribution 
is sensitive to the weights of orbits that have relatively large 
values of J z , with the consequence that it matters whether 
the centrifugal potential contains L z or \L Z \ +~f,J z , and we 
find materially better fits to the distributions of both V4, and 
v z when 7=1 rather than zero. 

Regardless of the value of 7, the AA predicts a value 
for (vr) 1 / 2 that exceeds the true value by an amount that 
grows with \z\, being ~ 3.4 per cent at \z\ = 1.5 kpc. The 
AA yields a value of <r^ that lies very close to the true value 
for \z\ < 1.3 kpc, but exceeds the true value by ~ 10 percent 
at 2 1 = 2 kpc because the true value declines with \z\ at 
\z\ > 1.3 kpc, whereas that obtained from the AA does not. 
With 7 = 1 the AA yields a value for (vz) 1 ^ 2 that lies within 
3 percent of the true value right up to \z\ = 2.3 kpc. 

The AA inevitably predicts that the velocity ellipsoid 
has two axes parallel to the plane, so we must turn to the 
full torus model to discover how the velocity ellipsoid tilts 
as one moves away from the plane. We find that its longest 
axis points quite close to the Galactic centre. This result 
emerges through averaging the quite disparate contributions 
of individual tori. Consequently it reflects the structure of 
the df as much as the gravitational potential. 

From a computational perspective, the AA is extremely 
convenient, both because it does not require specialised 
torus-generating code, and because it yields J from (x,v) 
rather than (x, v) from J. Consequently, a model's observ- 
ables can be obtained by integrating over velocity space, just 
as traditionally we have obtained the observables of models 
with DFs of the form f(E) and f(E, L). While McMillan & 
Binney (2008) showed that it is possible to determine J from 
(x,v), the procedure used is iterative and time-consuming, 
so for this paper observables were estimated by binning the 
particles of a discrete realisation obtained by Monte-Carlo 
sampling the df. Even this procedure is computationally 
expensive when enough samples are drawn to make Pois- 
son noise negligible, so it is very useful to be able to obtain 
good approximations to J(x,v) from the AA, and we antic- 
ipate that the AA will be widely used in the interpretation 
of observations of the Galaxy. 

Paper I and the present paper represent two small steps 
towards the kind of Galaxy modelling apparatus that should 
available before a preliminary Gaia Catalogue appears in the 
second half of this decade. The next big step is to carry the 
predictions of models into the space of observables - such 
as apparent magnitudes, parallaxes and proper motions - 
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and then to explore how tightly the df of stars can be con- 
strained by data of varying extent and precision. This step 
is crucial because distance uncertainties propagate from ob- 
servables such as magnitudes and proper motions to corre- 
lated errors in physical quantities such as stellar masses and 
velocities. We hope to report on this work soon. 



REFERENCES 

Aumer M. & Binney J.J., 2009, MNRAS, 397, 1286 
Bahcall J.N., Soneira R.M., 1984, ApJS, 55, 67 
Binney J., 2010, MNRAS, 401, 2318 (Paper I) 
Binney J., Spergel D.N., 1984, MNRAS, 206, 159 
Binney J., Tremaine S., 2008, "Galactic Dynamics", 

Princeton University Press, Princeton (BT08) 
Dehnen W., Binney J., 1998, MNRAS, 294, 429 
Dehnen W., Gerhard O.E., 1993, MNRAS, 261, 311 
Gebhardt, K. et al (15 authors), 2003, ApJ, 583, 92 
Gerhard O.E., Saha P., 1991, MNRAS, 251, 449 
Gilmore G., Reid N., 1983, MNRAS, 202, 1025 
Hfner R., Evans N.W., Dehnen W., Binney J., 2000, MN- 
RAS, 314, 433 
Kaasalainen M., 1995, MNRAS, 275, 162 
Kalnajs A., 1977, ApJ, 212, 637 

Krajnovic D., Cappellari M., Emsellem E., McDermid 

R.M., de Zeeuw P.T., 2005, MNRAS, 357, 1113 
Malin D.F., Carter D., 1980, Nature, 285, 643 
McMillan, P., Binney J., 2008, MNRAS, 390, 429 
Miyamoto M., Nagai R., 1975, PASJ, 27, 533 
Quinn P., 1984, ApJ, 279, 596 

Robin A.C., Reyle C, Derriere S., Picaud, S., 2003, A&A, 

409, 523 
Rowley G., 1988, ApJ, 331, 124 
Schonrich R., Binney J., 2009, MNRAS, 396, 203 
Shu F.H., 1969, ApJ, 158, 505 
Skrutskic M.F., et al., 2006, AJ, 131, 1163 
Schwarzschild M., 1979, ApJ, 232, 236 
Schweizer F., Seitzer P., 1992, AJ, 104, 1039 
Sharma, S., Bland-Hawthorn J., Johnston K.V., Binney J., 

2010, ApJ xxx 
Siebert, A. et al. (19 authors) 2008, MNRAS, 391, 793 
Smith M.C., Evans N.W., An J.H., 2009, ApJ, 698, 1110 
Thomas J., Saglia R.P., Bender R., Thomas D., Gebhardt 
K., Magorrian J., Corsini E.M., Wegner G., 2005, MN- 
RAS, 360, 1355 
Weinberg, M., 1994, ApJ, 421, 481 



APPENDIX A: EVALUATING JACOBIANS 

The observables of individual orbits can be obtained from 
the df /(J) = S(J — Jo), where Jo gives the orbits actions. 
Then 



(v i v j ) x = J d\v i v j f(3) = J d 3 J ggj 



ViVj8(J — Jo) 

(11) 



where 8i are the phases that bring the star to x. 

Since both the (x, v) and (8, J) systems are canonical, 



there exists a generating function G(J,x) for the transfor- 
mation (x, v) -o- {J, 9). So 

9 - — v - — 

dJi ' 3 dxj 

Differentiating again 

dOA 8 2 G 
dxj I dxjdJi 



(dVj\ 

\dJi), 



(12) 



(13) 



Taking determinants 



d(8) 



9(x) 



9(v) 



8(3) 



(14) 



The Jacobian on the left is the orbit's probability density 
in real space because the amount of time a star spends in 
a region of angle space is proportional to its volume, d 3 #. 
It makes good physical sense that the Jacobian in equation 
(11) is the orbit's probability density. 

With toy variables distinguished from true ones by a 
superscript T, the generating function of the transforma- 
tion (8, J) -H> (8^ T \ j( T )) between true and toy angle-action 
variables is 



S(J, (T) ) = J ■ 6> (T) + ^ Sn(J) sin(n ■ 8 m ), 
1 1 

where n is a vector with integer components, so 
J (T) = J + ^s n ncos(n - (T) ) 



(15) 



8 = 8 m +Y %sin(n.0 



3< T h 



53 



(16) 



The torus machine delivers the numerical values of both s 
and 8s n /83. We need 

(d(y)\ = ( 3(x) \ ( 8(8^) 
\d(8)J, \d(8^)JA 9(8) 



(17) 



J )/ j \ v ' / J 

The inverse of the second Jacobian on the right follows triv- 
ially from equation (16), and for the first Jacobian we can 
write 



dx = 



/ 8y 



(T) VaJ( T )/e< T ) 



(18) 



Dividing through by dd^ and holding J constant we find 
/ d(x) \ f dx \ f 9x \ ( 8J (T) \ 

\d{6W)) 3 ~ Ue (T V J( T) + Ujmjflcr, ' [dd^Jj 

The first two matrices on the right involve only toy variables 
so they are available analytically, and the third matrix can 
be obtained from equations (16). 



.(19) 



